Wilms tumor, SRP358696 bulk RNAseq analysis

Preliminary analysis of the SRP358696 dataset, which contain Wilms tumor samples.

set.seed(123)

library(DESeq2)
library(fgsea)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(EnhancedVolcano)
library(ggvenn)
source("/storage/research/dbmr_rubin_lab/scripts/save_png_pdf.R")
source("/storage/research/dbmr_rubin_lab/scripts/fgsea_save_res.R")

plotdir=("plots/")
if(!dir.exists(plotdir)){dir.create(plotdir,  recursive = T)}

Load sample names

samples <- scan("../samples.txt", character())

Load rsem counts

rsem_res_coding <- read.table("../rsem/all.genes.expected_count.results_coding", header = T, sep = "\t", row.names=1, check.names = F)

rsem_res_coding$gene_name_uniq <- make.unique(rsem_res_coding$gene_name)

remove metadata cols (note new rsem has more cols) last col is gene_name_uniq

counts <- rsem_res_coding[c(-1:-27, -ncol(rsem_res_coding))]

rownames(counts) <- rsem_res_coding$gene_name_uniq

counts need to be integers

counts <- round(counts)

remove low counts

counts <- counts[rowSums(counts) > 5,]

Load metadatra

meta <- openxlsx::read.xlsx("../GSE200256_metadata_bulk_RNA_seq.xlsx")

meta$sample <- sub("_", "", meta[["Sample.Name"]])

assure same order

identical(meta$sample, colnames(counts))

[1] TRUE

DESeq2

set coldata

coldata <- data.frame(condition = gsub(" .+", "", meta$isolate),
                        tissue = meta$tissue,
                        sex = meta$sex,
                    row.names =  colnames(counts))

construct a DESeqDataSet

dds <- DESeqDataSetFromMatrix(countData = counts,
                            colData = coldata,
                            design = ~ condition + tissue + sex)
#relevel
relevel(dds$condition, ref = "Favorable")

run DGE inference

dds <- DESeq(dds)

Variance stabilizing transformation

vsd <- vst(dds, blind=FALSE)

z-scores of vst data (for visualisation etc.)

vsd.zscore <- as.data.frame(t(scale(t(assay(vsd)), scale=TRUE, center=TRUE)))

Heatmap of the sample-to-sample distances

sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)

colors <- colorRampPalette(rev(brewer.pal(9, "Blues")) )(255)

# use fixed cell sizes to look symmetrical, then adjust figure size to have proper margins
p <- pheatmap::pheatmap(sampleDistMatrix,
        clustering_distance_rows=sampleDists,
        clustering_distance_cols=sampleDists,
        col=colors,
        annotation = coldata,
        main = "Sample distance\n",
        cellwidth = 12, cellheight = 12,
        silent = T)

save_png_pdf(p, paste0(plotdir, "heatmap"), height = 12, width = 12)

PCA

pca <- list("PC1-PC2" = 1:2,
            "PC3-PC4" = 3:4,
            "PC5-PC6" = 5:6,
            "PC7-PC8" = 7:8,
            "PC9-PC10" = 9:10)

for (x in names(pca)) {
    cat(x, "\n")
    pcaData <- plotPCA(vsd, intgroup=c("condition"), returnData=TRUE, pcsToUse = pca[[x]])
    percentVar <- round(100 * attr(pcaData, "percentVar"))
    
    p <- ggplot(pcaData, aes(!!sym(paste0("PC", pca[[x]][1])), !!sym(paste0("PC", pca[[x]][2])), fill=condition, color=tissue, shape=sex)) +
        geom_point(size=3, stroke = 1.2) +
        scale_shape_manual(values = 21:25) + # use fillable shapes
        scale_color_brewer(type = "qual", palette = "Set1",
                        guide = guide_legend(override.aes = list(shape = 1))) + # use a hollow shape in legend
        scale_fill_manual(values = c("Anaplastic" = "grey20", "Favorable" = "grey90"),
                        guide = guide_legend(override.aes = list(shape = 23, stroke=0.5))) + # use fillable shape in legend
        xlab(paste0("PC",pca[[x]][1],": ",percentVar[1],"% variance")) +
        ylab(paste0("PC",pca[[x]][2],": ",percentVar[2],"% variance")) +
        xlim(c(-max(abs(pcaData[[paste0("PC", pca[[x]][1])]]))-1, max(abs(pcaData[[paste0("PC", pca[[x]][1])]]))+1)) + # PC can be small that dots get trimmed. Expand a bit.
        ylim(c(-max(abs(pcaData[[paste0("PC", pca[[x]][2])]]))-1, max(abs(pcaData[[paste0("PC", pca[[x]][2])]]))+1)) + # PC can be small that dots get trimmed. Expand a bit.
        geom_text_repel(aes(label = name),
                        size = 1.5,
                        segment.size = 0.1,
                        colour = 'grey50',
                        box.padding   = 0.10,
                        segment.color = 'grey50',
                        force = 50,
                        max.overlaps = 30,
                        show.legend = F) +
        theme_minimal()
    
    save_png_pdf(p, paste0(plotdir, x), height = 5, width = 8)
    print(p)
}

PC1-PC2 
PC3-PC4 
PC5-PC6 
PC7-PC8 
PC9-PC10 

verify comparisons

resultsNames(dds)

[1] "Intercept"                         "condition_Favorable_vs_Anaplastic" "tissue_Kidney_vs_Diaphram"        
[4] "tissue_Left.Kidney_vs_Diaphram"    "tissue_Liver_vs_Diaphram"          "tissue_Right.Kidney_vs_Diaphram"  
[7] "sex_male_vs_female"               

Differential expression results for each comparison There is only 1 comparison

c <- "Anaplastic_vs_Favorable"
res <- results(dds, contrast = c("condition", "Anaplastic", "Favorable"), alpha = 0.05)

MA-plot

save_png_pdf(plotMA(res, ylim=c(-10,10)), paste0(plotdir, "DESeq2_res.", c, ".MA-plot"))

res to table, and add gene info

res.df <- as.data.frame(res)


res.df <- merge(res.df, rsem_res_coding[,c("gene_name_uniq", "gene_name", "gene_id", "gene_type", "seqnames", "start", "end", "width", "strand")], by.x=0, by.y="gene_name_uniq", all.x = T, all.y = F, sort = F)
colnames(res.df)[1] <- "gene_name_uniq"

save results

openxlsx::write.xlsx(res.df, file="DESeq2_res.xlsx")

volcano plots

EnhancedVolcano

p <- EnhancedVolcano(res.df,
                    lab = res.df$gene_name,
                    x = 'log2FoldChange',
                    y = 'padj',
                    ylab = bquote(~-Log[10] ~ italic(Padj)),
                    title = c,
                    subtitle = NULL,
                    cutoffLineCol = "gray10",
                    cutoffLineWidth = 0.2,
                    pointSize = 2,
                    labSize = 3,
                    axisLabSize = 15,
                    titleLabSize = 15,
                    labCol = "gray30",
                    col = c("grey60", "slategray3", "lightpink", "tomato"),
                    legendPosition = 'right',
                    legendLabels=c("NS", expression(Log[2] ~ FC), "Padj", expression(Padj ~ and ~ Log[2] ~ FC)),
                    legendLabSize = 10,
                    max.overlaps = 50,
                    drawConnectors = T,
                    widthConnectors = 0.1,
                    arrowheads = F,
                    colConnectors = "grey30",
                    directionConnectors="both" # 'y' works ok for fewer genes. ideal should be upregulated to one side, down to the other
                    ) +
    theme_minimal()

save_png_pdf(p, paste0(plotdir, "EnhancedVolcano.", c), height = 6, width = 8)

fgsea Msigdb

Load get gene sets in gmt

msigdb.hs.gmt <- readRDS("/storage/research/dbmr_rubin_lab/resources/msigdb/msigdb.hs.gmt.rds")
msigdb.hs.info <- readRDS("/storage/research/dbmr_rubin_lab/resources/msigdb/msigdb.hs.info.df.rds")

make rank

res.rank <- setNames(res.df$stat, make.unique(res.df$gene_name))

run only relevant collections

collections <- c("C2.CP", "C2.CP:KEGG_LEGACY", "C2.CP:KEGG_MEDICUS", "C2.CP:REACTOME", "C4.3CA", "C4.CGN", "C4.CM", "C5.GO:BP", "C5.GO:MF", "C6", "C7.IMMUNESIGDB", "C8", "H")

run fgsea

res.fgseaRes.msigdb <- sapply(collections, function(x) {
    cat(c, x, "\n")
    fgsea(pathways = msigdb.hs.gmt[[x]],
            stats = res.rank,
            eps = 0.0, minSize = 15, maxSize = 500)
}, simplify=F)

# Save session
save.image(file="session.RData")

Plot fgsea results

fgsea_save_res(res.fgseaRes.msigdb,
                rank = res.rank,
                basename = "fgsea.",
                gmt = msigdb.hs.gmt[collections],
                gmt_info = msigdb.hs.info,
                suffix = c,
                subtitle = c)

DE on transcript level

load rsem counts

rsem_res_coding.transcript <- read.table("../rsem/all.isoforms.expected_count.results_coding", header = T, sep = "\t", row.names=1, check.names = F)

# remove metadata cols (note new rsem has more cols)
counts.transcript <- rsem_res_coding.transcript[-1:-27]

# rownames can be transcript_name (they are unique)
sum(duplicated(rsem_res_coding.transcript$transcript_name))
# [1] 0
# but replace "-" with "_" so DESeq2 does not complain
rownames(counts.transcript) <- sub("-", "_", rsem_res_coding.transcript$transcript_name)

Prepare counts

# counts.transcript need to be integers
counts.transcript <- round(counts.transcript)

# remove low counts.transcript
counts.transcript <- counts.transcript[rowSums(counts.transcript) > 5,]

DESeq2

construct a DESeqDataSet

dds.transcript <- DESeqDataSetFromMatrix(countData = counts.transcript,
                            colData = coldata,
                            design = ~ condition + tissue + sex) 

# relevel
relevel(dds.transcript$condition, ref = "Favorable")

run DGE inference

dds.transcript <- DESeq(dds.transcript)

Variance stabilizing transformation

vsd.transcript <- vst(dds.transcript, blind=FALSE)

z-scores of vst data (for visualisation etc.)

vsd.transcript.zscore <- as.data.frame(t(scale(t(assay(vsd.transcript)), scale=TRUE, center=TRUE)))

Heatmap of the sample-to-sample distances

sampleDists <- dist(t(assay(vsd.transcript)))
sampleDistMatrix <- as.matrix(sampleDists)

colors <- colorRampPalette(rev(brewer.pal(9, "Blues")) )(255)

# use fixed cell sizes to look symmetrical, then adjust figure size to have proper margins
p <- pheatmap::pheatmap(sampleDistMatrix,
        clustering_distance_rows=sampleDists,
        clustering_distance_cols=sampleDists,
        col=colors,
        annotation = coldata,
        main = "Sample distance\n",
        cellwidth = 12, cellheight = 12,
        silent=T)

save_png_pdf(p, paste0(plotdir, "heatmap.transcript"), height = 12, width = 12)

PCA

pcaData <- plotPCA(vsd.transcript, intgroup=c("condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

p <- ggplot(pcaData, aes(PC1, PC2, fill=condition, color=tissue, shape=sex)) +
    geom_point(size=3, stroke = 1.2) +
    scale_shape_manual(values = 21:25) + # use fillable shapes
    scale_color_brewer(type = "qual", palette = "Set1",
                    guide = guide_legend(override.aes = list(shape = 1))) + # use a hollow shape in legend
    scale_fill_manual(values = c("Anaplastic" = "grey20", "Favorable" = "grey90"),
                    guide = guide_legend(override.aes = list(shape = 23, stroke=0.5))) + # use fillable shape in legend
    xlab(paste0("PC1: ",percentVar[1],"% variance")) +
    ylab(paste0("PC2: ",percentVar[2],"% variance")) +
    ylim(c(-max(abs(pcaData$PC2))-1, max(abs(pcaData$PC2))+1)) + # if PC2 is small dots get trimmed. Expand a bit.
    geom_text_repel(aes(label = name),
                    size = 1.5,
                    segment.size = 0.1,
                    colour = 'grey50',
                    box.padding   = 0.10,
                    segment.color = 'grey50',
                    force = 50,
                    max.overlaps = 30,
                    show.legend = F) +
    theme_minimal()

save_png_pdf(p, paste0(plotdir, "PCA.transcript"), height = 5, width = 8)

null device 
          1 

print(p)

Fetch results

res.transcript <- results(dds.transcript, contrast = c("condition", "Anaplastic", "Favorable"), alpha = 0.05)

MA-plot

save_png_pdf(plotMA(res.transcript, ylim=c(-20,20)), paste0(plotdir, "DESeq2_res.transcript.", c, ".MA-plot"))

res to table, and add gene info

res.transcript.df <- as.data.frame(res.transcript)

# _ to - so that we can merge
res.transcript.df$transcript_name <- sub("_", "-", rownames(res.transcript))

res.transcript.df <- merge(res.transcript.df, rsem_res_coding.transcript[,c("transcript_name", "transcript_id", "transcript_type", "gene_name", "gene_id", "seqnames", "start", "end", "width", "strand")], by="transcript_name", all.x = T, all.y = F, sort = F)

save results

openxlsx::write.xlsx(res.transcript.df, file="DESeq2_res.transcript.xlsx")

volcano plots

EnhancedVolcano

p <- EnhancedVolcano(res.transcript.df,
                    lab = res.transcript.df$transcript_name,
                    x = 'log2FoldChange',
                    y = 'padj',
                    ylab = bquote(~-Log[10] ~ italic(Padj)),
                    title = c,
                    subtitle = NULL,
                    cutoffLineCol = "gray10",
                    cutoffLineWidth = 0.2,
                    pointSize = 2,
                    labSize = 2,
                    axisLabSize = 15,
                    titleLabSize = 15,
                    labCol = "gray30",
                    col = c("grey60", "slategray3", "lightpink", "tomato"),
                    legendPosition = 'right',
                    legendLabels=c("NS", expression(Log[2] ~ FC), "Padj", expression(Padj ~ and ~ Log[2] ~ FC)),
                    legendLabSize = 10,
                    max.overlaps = 50,
                    drawConnectors = T,
                    widthConnectors = 0.1,
                    arrowheads = F,
                    colConnectors = "grey30",
                    ) +
    theme_minimal()

save_png_pdf(p, paste0(plotdir, "EnhancedVolcano.transcript.", c), height = 6, width = 8)

Kallisto

load kallisto counts over transcripts

kallisto_res.transcript <- read.table("../kallisto/all.est_count.txt.gz", header = T, sep = "\t", row.names=1, check.names = F)[-1] # omit column "length"

use tximport to aggregate to gene level

Prepare data

read the gtf used for the Kallisto index (full anno, not just “annotation” or “primary_assembly.annotation”)

gencode.v48 <- rtracklayer::readGFF("/storage/research/dbmr_rubin_lab/pipeline/ref/anno/hg38/gencode.v48.chr_patch_hapl_scaff.annotation.gtf.gz")

did we fetch everything?

sum(!rownames(kallisto_res.transcript) %in% gencode.v48$transcript_id)

[1] 0

more sanity checks

rsem_transcript_list <- read.table("../rsem/all.isoforms.expected_count.results", header = T, sep = "\t")[[1]]
sum(rsem_transcript_list %in% gencode.v48$transcript_id)

[1] 387954

length(rsem_transcript_list)

[1] 387954

make a data.frame called tx2gene with two columns: 1) transcript ID and 2) gene ID. The column names do not matter but this column order must be used.

tx2gene <- gencode.v48[c("transcript_id", "gene_id")]

# remove NA and duplicate transcripts
tx2gene <- na.omit(tx2gene)
tx2gene <- tx2gene[!duplicated(tx2gene$transcript_id),]

gene counts from kallisto TSV files

k_files <- file.path("../kallisto", samples, "abundance.tsv.gz")
names(k_files) <- samples

tximport to aggregate to gene level

# we don't need ignoreAfterBar
kallisto_res.gene <- tximport::tximport(k_files, type = "kallisto", tx2gene = tx2gene)

# to df
kallisto_res.gene.df <- as.data.frame(kallisto_res.gene[["abundance"]])

compare rsem and kallisto

rsem_res <- read.table("../rsem/all.genes.expected_count.results", header = T, sep = "\t", row.names=1, check.names = F)[c(-1:-27)] # omit metadata cols

rsem_res.transcript <- read.table("../rsem/all.isoforms.expected_count.results", header = T, sep = "\t", row.names=1, check.names = F)[c(-1:-27)] # omit metadata cols

sanity check

identical(colnames(rsem_res.transcript), colnames(kallisto_res.transcript))

[1] TRUE

Note: kallisto has more genes/transcripts than rsem (which are all in kallisto)

rsem_vs_kallisto <- list(transcript = cor(rsem_res.transcript, kallisto_res.transcript[rownames(rsem_res.transcript),], method = "spearman"))

rsem_vs_kallisto[["transcript.coding"]] <- cor(rsem_res_coding.transcript[-1:-27], kallisto_res.transcript[rownames(rsem_res_coding.transcript),], method = "spearman")

rsem_vs_kallisto[["gene"]] <- cor(rsem_res, kallisto_res.gene.df[rownames(rsem_res),], method = "spearman")

rsem_vs_kallisto[["gene.coding"]] <- cor(rsem_res[rownames(rsem_res_coding),], kallisto_res.gene.df[rownames(rsem_res_coding),], method = "spearman")

Plot correlatoions

for (x in names(rsem_vs_kallisto)) {
    p <- ComplexHeatmap::pheatmap(rsem_vs_kallisto[[x]],
        cluster_rows = F,
        cluster_cols = F,
        main = paste0("RSEM vs. Kallisto", "\n", "counts", "\n", x),
        fontsize_row = 8,
        fontsize_col = 8,
        cellheight = 12,
        cellwidth = 12,
        display_numbers = T,
        fontsize_number = 4,
        heatmap_legend_param = list(title = "Spearman"))

    save_png_pdf(p, paste0(plotdir, "rsem_vs_kallisto.", x), height = 10, width = 10)
    print(p)
}

Save session

save.image(file="session.RData")

sessioninfo::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.5.1 (2025-06-13)
 os       Ubuntu 24.04.3 LTS
 system   x86_64, linux-gnu
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Etc/UTC
 date     2026-02-06
 rstudio  2025.09.2+418 Cucumberleaf Sunflower (server)
 pandoc   3.8.2.1 @ /usr/bin/ (via rmarkdown)
 quarto   1.7.32 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package              * version  date (UTC) lib source
 abind                  1.4-8    2024-09-12 [1] RSPM (R 4.5.0)
 Biobase              * 2.68.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 BiocGenerics         * 0.54.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 BiocParallel           1.42.1   2025-06-01 [1] Bioconductor 3.21 (R 4.5.1)
 bslib                  0.9.0    2025-01-30 [2] RSPM (R 4.5.0)
 cachem                 1.1.0    2024-05-16 [2] RSPM (R 4.5.0)
 circlize               0.4.16   2024-02-20 [1] RSPM (R 4.5.0)
 cli                    3.6.5    2025-04-23 [2] RSPM (R 4.5.0)
 clue                   0.3-66   2024-11-13 [1] RSPM (R 4.5.0)
 cluster                2.1.8.1  2025-03-12 [3] CRAN (R 4.5.1)
 codetools              0.2-20   2024-03-31 [3] CRAN (R 4.5.1)
 colorspace             2.1-1    2024-07-26 [1] RSPM (R 4.5.0)
 ComplexHeatmap         2.24.1   2025-06-25 [1] Bioconductor 3.21 (R 4.5.1)
 cowplot                1.2.0    2025-07-07 [1] RSPM (R 4.5.0)
 crayon                 1.5.3    2024-06-20 [2] RSPM (R 4.5.0)
 data.table             1.17.8   2025-07-10 [2] RSPM (R 4.5.0)
 DelayedArray           0.34.1   2025-04-17 [1] Bioconductor 3.21 (R 4.5.1)
 DESeq2               * 1.48.2   2025-08-27 [1] Bioconductor 3.21 (R 4.5.1)
 digest                 0.6.37   2024-08-19 [2] RSPM (R 4.5.0)
 doParallel             1.0.17   2022-02-07 [1] RSPM (R 4.5.0)
 dplyr                  1.1.4    2023-11-17 [2] RSPM (R 4.5.0)
 EnhancedVolcano      * 1.26.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 evaluate               1.0.5    2025-08-27 [2] RSPM (R 4.5.0)
 farver                 2.1.2    2024-05-13 [2] RSPM (R 4.5.0)
 fastmap                1.2.0    2024-05-15 [2] RSPM (R 4.5.0)
 fastmatch              1.1-6    2024-12-23 [1] RSPM (R 4.5.0)
 fgsea                * 1.34.2   2025-07-13 [1] Bioconductor 3.21 (R 4.5.1)
 foreach                1.5.2    2022-02-02 [1] RSPM (R 4.5.0)
 fs                     1.6.6    2025-04-12 [2] RSPM (R 4.5.0)
 generics             * 0.1.4    2025-05-09 [2] RSPM (R 4.5.0)
 GenomeInfoDb         * 1.44.1   2025-07-23 [1] Bioconductor 3.21 (R 4.5.1)
 GenomeInfoDbData       1.2.14   2025-08-18 [1] Bioconductor
 GenomicRanges        * 1.60.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 GetoptLong             1.0.5    2020-12-15 [1] RSPM (R 4.5.0)
 ggplot2              * 4.0.0    2025-09-11 [2] RSPM (R 4.5.0)
 ggplotify            * 0.1.2    2023-08-09 [1] RSPM (R 4.5.0)
 ggrepel              * 0.9.6    2024-09-07 [1] RSPM (R 4.5.0)
 ggvenn               * 0.1.19   2025-10-08 [1] RSPM
 GlobalOptions          0.1.2    2020-06-10 [1] RSPM (R 4.5.0)
 glue                   1.8.0    2024-09-30 [2] RSPM (R 4.5.0)
 gridGraphics           0.5-1    2020-12-13 [1] RSPM (R 4.5.0)
 gtable                 0.3.6    2024-10-25 [2] RSPM (R 4.5.0)
 htmltools              0.5.8.1  2024-04-04 [2] RSPM (R 4.5.0)
 httr                   1.4.7    2023-08-15 [2] RSPM (R 4.5.0)
 IRanges              * 2.42.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 iterators              1.0.14   2022-02-05 [1] RSPM (R 4.5.0)
 jquerylib              0.1.4    2021-04-26 [2] RSPM (R 4.5.0)
 jsonlite               2.0.0    2025-03-27 [2] RSPM (R 4.5.0)
 knitr                  1.50     2025-03-16 [2] RSPM (R 4.5.0)
 labeling               0.4.3    2023-08-29 [2] RSPM (R 4.5.0)
 lattice                0.22-7   2025-04-02 [3] CRAN (R 4.5.1)
 lifecycle              1.0.4    2023-11-07 [2] RSPM (R 4.5.0)
 locfit                 1.5-9.12 2025-03-05 [1] RSPM (R 4.5.0)
 magrittr               2.0.4    2025-09-12 [2] RSPM (R 4.5.0)
 Matrix                 1.7-3    2025-03-11 [3] CRAN (R 4.5.1)
 MatrixGenerics       * 1.20.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 matrixStats          * 1.5.0    2025-01-07 [1] RSPM (R 4.5.0)
 patchwork            * 1.3.2    2025-08-25 [1] RSPM (R 4.5.0)
 pheatmap               1.0.13   2025-06-05 [1] RSPM (R 4.5.0)
 pillar                 1.11.1   2025-09-17 [2] RSPM (R 4.5.0)
 pkgconfig              2.0.3    2019-09-22 [2] RSPM (R 4.5.0)
 png                    0.1-8    2022-11-29 [2] RSPM (R 4.5.0)
 R6                     2.6.1    2025-02-15 [2] RSPM (R 4.5.0)
 rappdirs               0.3.3    2021-01-31 [2] RSPM (R 4.5.0)
 RColorBrewer         * 1.1-3    2022-04-03 [2] RSPM (R 4.5.0)
 Rcpp                   1.1.0    2025-07-02 [2] RSPM (R 4.5.0)
 rjson                  0.2.23   2024-09-16 [1] RSPM (R 4.5.0)
 rlang                  1.1.6    2025-04-11 [2] RSPM (R 4.5.0)
 rmarkdown              2.30     2025-09-28 [2] RSPM (R 4.5.0)
 rsconnect              1.6.0    2025-10-28 [2] RSPM (R 4.5.0)
 rstudioapi             0.17.1   2024-10-22 [2] RSPM (R 4.5.0)
 S4Arrays               1.8.1    2025-06-01 [1] Bioconductor 3.21 (R 4.5.1)
 S4Vectors            * 0.46.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 S7                     0.2.0    2024-11-07 [1] RSPM (R 4.5.0)
 sass                   0.4.10   2025-04-11 [2] RSPM (R 4.5.0)
 scales                 1.4.0    2025-04-24 [2] RSPM (R 4.5.0)
 sessioninfo            1.2.3    2025-02-05 [2] RSPM (R 4.5.0)
 shape                  1.4.6.1  2024-02-23 [1] RSPM (R 4.5.0)
 SparseArray            1.8.1    2025-07-23 [1] Bioconductor 3.21 (R 4.5.1)
 SummarizedExperiment * 1.38.1   2025-04-30 [1] Bioconductor 3.21 (R 4.5.1)
 tibble                 3.3.0    2025-06-08 [2] RSPM (R 4.5.0)
 tidyselect             1.2.1    2024-03-11 [2] RSPM (R 4.5.0)
 UCSC.utils             1.4.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 vctrs                  0.6.5    2023-12-01 [2] RSPM (R 4.5.0)
 withr                  3.0.2    2024-10-28 [2] RSPM (R 4.5.0)
 xfun                   0.53     2025-08-19 [2] RSPM (R 4.5.0)
 XVector                0.48.0   2025-04-15 [1] Bioconductor 3.21 (R 4.5.1)
 yaml                   2.3.10   2024-07-26 [2] RSPM (R 4.5.0)
 yulab.utils            0.2.1    2025-08-19 [1] RSPM (R 4.5.0)

 [1] /storage/research/dbmr_rubin_lab/R_lib/Rocker_rstudio_4.5.1
 [2] /usr/local/lib/R/site-library
 [3] /usr/local/lib/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────